AD-A139  315 
Unclassified 


ESTIMATING  TIME  AVERAGES  VIA  RANDOMLY-SPACED  1/1 

OBSERVAT IONS(U)  WISCONSIN  UNIV-MADISON  MATHEMATICS 
RESEARCH  CENTER  B  L  FOX  ET  AL .  JAN  84  MRC-TSR-2638 
DAAG29-80-C- 004 1  F/G  12/1  NL 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  Of  STANDARDS  1963-A 


Mathematics  Research  Center 
University  of  Wisconsin— Madison 
610  Walnut  Street 
Madison,  Wisconsin  53705 


January  1984 


(Received  October  14,  1983) 


DT1C 

SELECTE 

MAR  2  21984 

B 


Approved  for  public  release 
Distribution  unlimited 


Sponsored  by 

U.  S.  Army  Research  Office 
P.  0.  Box  12211 
Research  Triangle  Park 
North  Carolina  27709 

OTIC  FILE  COPY 


■°i  o 


Oi 


A 


84 


UNIVERSITY  OF  WI8C0N8IN-MADI80N 
MATHEMATICS  RESEARCH  CENTER 


ESTIMATING  TIME  AVERAGES  VIA  RANDOMLY-SPACED  OBSERVATIONS 
Bennett  L.  Fox*'1  and  Peter  W.  Glynn ^ 

\ 

Technical  Summary  Report  #2638 
January  1984 

ABSTRACT 

To  estimate  continuous-time  averages  via  randomly-spaced  observations  of 

^  A  *  ,  -'v-  '  1 

discrete-event  systems,  develop  a  point-process  framework  and  use  it  to 

generalize  both  regenerative  and  stationary-process  oriented  simulation 

r -  ? 

methodologies.  give  consistent  estimators,  central  limit  theorems,  and  an 
effective  bias-reducing  jackknife.  The  impact  on  indirect  estimation  of 
transaction  (customer)  averages  is  discussed. 


AMS  (MOS)  Subject  Classifications:  62M15,  65C05,  68J05,  60G55 
Key  Words:  Simulation,  point  processes 

Work  Unit  Number  5  (Optimization  and  Large  Scale  Systems) 


♦Dipartement  d’informatique  et  de  recherche  opirationnelle,  University  de 
Montrial,  Montrial,  Quibec  H3C  3J7 


Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041. 
Research  also  supported  by 

'The  Natural  Sciences  and  Engineering  Research  Council  of  Canada. 

2The  Graduate  School  of  the  University  of  Wisconsin-Madison. 


o 


SIGNIFICANCE  AND  EXPLANATION 

vln  many  stochastic  systems,  one  is  interested  in  estimating  steady-state 
expected  values.  When  Monte  Carlo  simulation  is  used  to  estimate  such 
parameters,  an  assessment  of  accuracy,  in  the  form  of  confidence  intervals,  is 
often  required.  Most  procedures  for  producing  such  confidence  intervals 
require  that  the  simulation  be  sampled  so  that  the  time  increments  between 
observations  are  all  equal.  This  is  difficult  to  accon^lish  in  a  discrete- 
event  simulation,  since  the  clock  which  drives  the  simulation  is  incremented 

in  a  random  fashion. ^Our  purpose,  in  this  paper,  is  to  show  how  methods  for 

\ 

dealing  with  equally  spaced  observations  can  be  adapted  to  run  on  the  random 
time  scale  of  the  driving  clock  for  the  simulation. 
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ESTIMATING  TIME  AVERAGES  VIA  RANDOMLY-SPACED  OBSERVATIONS 


*  1  2 
Bennett  L.  Fox  '  and  Peter  W.  Glynn* 

1,  INTRODUCTION . 

Let  0  •  Tq  <  Tj  <  ...  be  event  tinea ■  Though  the  associated  sequence 
with  Xy  e  R  is  not  necessarily  an  imbedded  Markov  chain,  we  call  X^  the  state  at 
time  Ty  -  somewhat  abusing  the  term.  To  define  the  state  X(t)  at  an  arbitrary  time 
t,  interpolate! 

fl» 

(1.1)  X(t>  -  l  Zl  (t> 

k  0  Tk'Tk+1 ' 

where  the  indicator  IA  is  1  or  0  depending  on  whether  or  not  t  e  A.  For  this 
definition  to  make  sense,  every  state  change  must  correspond  to  an  event  time.  The  state 
does  not  change  continuously.  It  jumps  at  discrete  (possibly  random)  times.  In  other 
words,  we  have  a  discrete-event  system.  Let  f  be  a  real-valued  function.  Put 

(1.2)  r(t)  -  J  f(X(s))ds  . 

c  0 

We  solve  the  steady-state  problem  estimate  the  limit  (when  it  exists) 

(1.3)  r  -  lim  r(t> 

t*w 

and  construct  confidence  intervals  for  r. 

To  do  this,  we  develop  a  point-process  framework  and  use  it  to  generalise  both 
regenerative  and  stationary-process  oriented  simulation  methodologies.  Simply  averaging 
the  ^y’s  generally  inconsistently  estimates  r.  The  Ty’s  are  not  necessarily 
regeneration  times.  We  work  with  generally  dependent  observations,  in  contrast  to 
regenerative  approaches. 
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1.1  Shifts.  We  use  the  set  of  functions  u  :  I0,*»)  ♦  Rd  that  are  right  continuous  and 
have  left  limits  to  describe  the  sample  space  Q.  Define  X  via 

(1.4)  X(*)  «  X(*,u)  -  !!»<•) 

for  u>  e  n.  For  any  random  variable  R  :  fl  ♦  (0,«),  define  a  right  shift  via 

(1.5)  eR  »  0R(M)(“)  “  «(R(»)  +  •)  • 

Let  A  c  *d,  put  s0  «  o,  and  set 

(1.6)  Sy+1  -  inf  ft  >  S*  :  X(t-)  *  X(t),  X(t)  e  A) 

where  X(t-)  is  the  left  limit  of  X  at  t.  If  A  "  Rd,  then  -  T^.  Define 

(1.7)  Pfx  e  •)  -  PfX  •  8e  e  •} 

n  s 

n 

and 

(1.8)  pfc{x  e  •}  -  p{x  •  9t  e  •} 

For  most  simulations  where  the  steady-state  limit  r  exists,  limit  probabilities  P  and 
P  exist  such  that  weak  convergence  holds: 

(1.9)  Pn  — >  P 
and 

(1.10)  Pt  — >  P 
where 

(1.11)  PfX  •  8  e  •>  -  PfX  E  •} 

S1 

and  (shift  invariance) 

(1.12)  PfX  see*}-  PfX  e  •} 
for  t  >  0. 

Example  1.1.  Let  X  be  a  delayed  (resp..  non-delayed)  regenerative  process  under 
*  * 

P  (resp.,  P).  Assume  that  X  regenerates  when  it  hits  A.  Then  Pn  ■  P  for  n  >  1 

though  generally  Pq  *  P> 

Example  1.2.  Let  X  be  as  in  example  1.1.  If  the  regenerative  process  is  positive 

recurrent  and  the  regeneration-spacing  distribution  has  a  non-trivial  Lebesgue  component, 

^  <* 

then  (1.10)  and  (1.12)  hold  (see  Miller  (1972)).  Also,  P  and  P  are  related  by 


-2- 


(  1. 13 ) 


p{x  £  •}  -  zr-  I  Ptt  *  8  e  •»  b  >  t}<Jt 

■  M  A  t  I 


ES1  0 


where  e  denotes  expectation  under  P.  In  Section  3  we  show  that  (1.13)  holds  store 
generally. 

A  process  X  satisfying 

(1.14)  p{x  •  e_  c  •)  -  p{x  e  •) 

si 

is  synchronous  with  respect  to  the  imbedded  point  process  sequence  (S^) . 

Example  1.3.  Suppose  that  {Xn>  is  a  stationary  sequence  with  Xn  “Tn  ■  n.  Then 
X(t)  is  synchronous  with  respect  to  {s^} . 

From  examples  1 . 1  and  1 . 3  we  see  that  synchronous  processes  generalise  both  non- 
delayed  regenerative  processes  and  stationary  sequences.  Assume  that  the  simulator 
(somehow)  chooses  the  time  origin  so  that  everything  representing  the  initial  "transient* 
phase  is  to  its  left.  This  is  trivial  for  regenerative  processes  but  not  in  general.  This 
deletion  assumption  translates  mathematically  asi  (1.14)  holds. 

Let 

(1.15)  °„  “  sn+1  ~  sn 
and 


(1.16) 


X  (t) 
n 


'  x(s  +  t)i  0  <  t  <  a 

n  n 

•  i  t  >  a  . 

n 


We  assume  that 
(1.17) 


and 

(1.18) 

where 


(X  }  is  a  ^-mixing  sequence  (Billingsley  (I960), 
n 

pp.  166*168)  with  +-*dxing  coefficients  satisfying 


<  • 


E(Y  (|f|)2  *  a)  <  m 
n  A 
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(1.19) 


Yn(f)  «  j  f(Xn(t))dt 
Sn+1 

-  J  f (X(t) )dt  . 

Sn 

The  cycle  sequence  {  X^}  is  not  necessarily  iid.  If  it  were,  we  could  simply  apply 
regenerative  methodology.  Three  more  definitions  set  the  stage: 


n-1 

(1.20)  5(f)  *  1  Y.ffJ/n 

k»0 

(1.21)  a  «  S  /n 

n  n 

(1.22)  r  -  Y  (f)/5  . 

n  n  n 

1.2,  Preview  of  results.  In  Section  2  we  assume  that  (1.14),  (1.17),  and  (1.18)  hold, 
usually  without  further  explicit  mention.  We  show  that 

(1.23)  rn  ♦  r 

(1.24)  r(t)  ♦  r 

(1.25)  r  -  EY0(f)/Ea0 

and  generalize  the  "inspection  paradox" 


(1.26) 


I[x,-)(aN(s)  )dS 


P{aQ  >  s}ds/EaQ  . 


The  left  side  of  (1.26)  is  the  proportion  of  time  over  ( 0 , t)  that  the  cycle  in  progress 
had  length  at  least  x.  Its  limit,  the  right  side  of  (1.26),  is  precisely  that  for  the 
regenerative  case:  e.g.,  see  Bratley,  Fox,  and  Schrage  [(1983),  problem  3.7.4)  or  Heyman 
and  Sobel  [(1982),  $5.5).  Thus,  for  synchronous  processes  the  cycle  in  progress  at  time 
t  tends  to  be  longer  than  average,  thereby  biasing  r(t). 

Next,  Section  2  proves  a  central  limit  theorem  (CLT) ; 

(1.27)  n^  <rn  -  r)  =~>  ON(0,1)/Ea0 

(1.28)  t/2  (r(t)  -  r)  -*>  <JN(0,1)/(Ect0)1/2 

-4- 


where  N(Oa1)  is  a  aero-mean,  unit-variance,  normal  random  variable, 
(1.29)  O2  -  BY  (f)2  +  2  l  EY.  (f  )V  (£ ) 


(1.30) 


f<x)  -  f(x) 


In  other  words,  Y^  ( f )  ■  Y^tf)  -  a^r.  If  the  synchronous  process  X  is  either 
regenerative  or  stationary  with  -  n,  the  CLT  simplifies.  For  the  former,  the  second 
term  in  (1.29)  -  corresponding  to  covariances  -  vanishes.  For  the  latter,  an  =  1.  Also 
(1.31)  when  A  -  ***, 

d2  -  Ef2<X0)r02  ♦  2  I  «V«VVk 

k-1 


Tk  “  Tk+1  "  Tk  • 

2 

To  estimate  a  in  the  regenerative  case,  one  uses  the  knowledge  that  the  covariance 

2 

terms  in  (1.29)  vanish,  to  construct  the  estimator  at 

n 

*1 "  l 

k-1 

2  2 
the  estimator  <»n  is  easily  shown  to  be  strongly  consistent  for  c  . 

2  2 

In  the  general  case,  estimation  of  o  is  more  complicated.  The  parameter  a  can 
be  expressed  in  the  form 


o2  -  [var  Y„(f)  +  2  \  cov(YQ (f ) .Y^tf ) ) ] 


r[cov(Y„(f ) ,a„)  +  2  cov^tf),^)] 

k-1 


-  r[cov(a0,Y0(f ) )  +  2  \  cov(aQ,Y  (f))] 


+  r  [var  a.  +  2  \  cov(a.,o  )] 

k-1  u  K 


-5- 


the  four  bracketed  terms  appearing  in  (1.32)  can  be  consistently  estimated  by 

c1n'c2n'c3n'c4n  (say)>  this  can  be  accomplished  by  standard  techniques  (e.g.  batch  means, 

spectral  methods,  autoregressive  procedures).  Such  methods  include  parameters  (e  g.  batch 

sire,  spectral  window,  autoregressive  order)  which  must  be  keyed  to  sang>le  sire.  In  any 
2 

case,  the  estimator  o  is  then  given  by 
n 

2  2 
o  “  c,_  -  r  c,  -  r  c,_  +  r„c.  . 
n  »n  n  2n  n  3n  n  4n 

When  using  r(t),  replace  n  by 

(1.33)  N(t)  =  #  Sj^’s  observed  in  (0,t) 

Glynn  and  Iglehart  (1981),  Jow  ((1982),  pp.  54-56],  and  Streller  ((1980),  Theorem  3.2] 

prove  essentially  (1.27)  but  under  significantly  different  hypotheses. 

From  (1.27)  and  (1.28)  we  get  the  respective  confidence  intervals 

[r  -  z  o  /a  /n,  r  +  r.o  /a  /n] 
n  fin  n  n  6nn 

and 

tr(t)  -  V»(t)/(tVt))l/2'  r(t)  +  *6°»(t)/(t“N(t),V2] 

where  the  percentile  z  is  the  unique  solution  of  P[N(0,1)  <  r  ]  «  1  -  S/2.  If  o  is 
S  S  n 

a  Strongly  consistent  estimator  for  a,  then  these  are  asymptotically  exact  100(1  -  6)% 
confidence  intervals  for  r. 

Under  an  additional  assumption 

(1.34  there  exists  K  >  0  such  that  P{a  <  K)  ■>  1 

0 

and  sup{|f(x)|  :  x  e  R  }  <  K 

Section  2  proves  that 

(1-35)  Ern  -  r  -  B/(Ea0>2n  +  0(1/n) 

and  finds  an  expression  for  6.  The  form  of  (1.35)  motivates  a  jackknife: 

(1-36)  r  =  2r2n  -  (r(0,n  -  1)  +  r(n,2n  -  1))/2 

where 

b  b 

(1.37)  r(a,b)  -  (  £  Y  (f ) )/(  J  a  ) 

j“«  3  j-a  3 

and 
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(1.38) 


r  +  0(  t/n) 


reducing  bias  by  an  order  of  magnitude. 

Sometimes  there  is  no  bias.  If  Sn  -  Tn  ■  n,  then  Ern  -  r. 

For  another  case,  let  A  *  and  suppose  that  the  T^'s  see  time  averages.  For  the 
latter,  suppose  that 

(i)  N(t)  is  a  stationary  Poisson  process 
(ii)  for  each  t  >  0,  {X(s-)  :  0  <  s  <  t}  is  independent  of 

{N(t  +  u)  -  N(t)  :  u  >  0}. 

Wolff  (1982)  shows  that,  if  r(t)  *  r  a.s.,  then 

-  i  n_1 

(1.39)  r  -  -  l  f(X  >  +  r  a.s. 

k“0 

A 

Since  the  summands  are  identically  distributed,  Er^  ”  r.  Section  2  concludes  by  proving 
that 

(1.40)  /2n  (r  -  r)  —>  OM(0,1)/Ea0  . 

Combining  (1.27)  and  (1.40),  we  see  that  rJn  and  r2n  have  the  same  asymptotic 
variance.  So  our  jackknife  reduces  bias  without  increasing  variance. 

Section  3  begins  by  proving  that 

,  t  - 

(1.41)  -  I  P{x  •  8  e  »}ds  ♦  P{X  e  •} 

t  o  8 

where  P  satisfies  (1.13),  still  assuming  that  (1.14),  (1.17),  and  (1.18)  hold.  The  point 

is  that  (1.13)  remains  valid  when  X  is  (merely)  synchronous.  Franken  et  al.  [(1982), 

Chapter  1]  show  that  it  holds  under  even  weaker  conditions. 

When  X  is  regenerative,  P{X  c  •}  “  P{X  •  8  e  •),  "inverting"  (1.13).  This 

S1 

simple  inversion  generally  fails,  assuming  only  ( 1 . 9 )  —  ( 1 . 12  ) .  Recall  that  under  P,  the 
cycle  trapping  a  fixed  time  tends  to  be  longer  than  a  typical  cycle.  In  the  regenerative 
case,  only  the  actual  trapping  cycle  is  affected,  due  to  independence  of  cycles.  In 
general,  however,  neighboring  cycles  are  also  affected  due  to  cycle  correlation.  Franken 
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et  al.  ((1982),  p.  23)  show  that  the  general  Inverse  to  (1.13)  is 

(1.42)  P(X  E  •)  =■  y  l  ?{X  •  9  e  •;  S.  <  1) 

A  k-1  sk  k 

where  3  “  EN(1)  and  E  denotes  expectation  under  P.  Comparing  (1.11)  and  (1.14),  we 
see  that  X  is  synchronous  under  P,  the  Palm  distribution  of  P.  Here  A  is  the 
intensity  of  P,  and  X  is  stationary  under  P. 

Chapter  1  of  Franken  et  al.  (1982)  thoroughly  discusses  the  relationships  between  P 
and  P  and  proves  an  intuitive  alternative  to  (1.42): 

(1.43)  P{X  e  •}  =*  lim  P{X  •  9  e  «|S  <  h)  . 

hiO  S1  1 

This  shows  that  the  Palm  distribution  P  is  the  stationary  diet-  cl  ->n  P  conditioned  on 
hitting  A  at  time  0. 

We  want  (1.12)  to  hold  with  P  replacing  P.  This  means  '  the  simulator  (somehow) 
deletes  the  entire  transient  phase,  choosing  the  time  origin  so  -la  nis  phase  is  to  itB 
left.  Typically  this  is  impossible  to  do  exactly  in  practice,  but  to  proceed 
mathematically  we  assuaie  it  has  been  done  exactly.  This  translates  as 

(1.44)  P{X  *  9t  e  •}  =  P{X  e  •) 

holds.  This  is  stronger  than  (1.14)  because  Er(t)  *  r  under  (1.44)  but  not  generally 
under  (1.14).  In  fact  (1.44)  usually  holds  literally  only  if  the  initial  state  is 
generated  by  the  (generally  unknown)  stationary  distribution  P.  Even  for  regenerative 
processes,  where  synchronization  via  (1.14)  is  trivial,  making  (1.44)  hold  even  to  a  first 
approximation  is  generally  hard.  Nevertheless,  in  the  rest  of  this  section  and  in  Section 
3,  usually  without  further  explicit  mention,  we  assume  that  (1.44)  holds  and  that,  under 
P  but  not  necessarily  under  P,  (1.14),  (1.1’),  and  (1.18)  hold.  By  contrast.  Section  2 
assumes  that  (1.14),  (1.17),  and  (1.18)  hold  under  P. 

Several  counterparts  to  results  in  Section  2  are  proved  in  Section  3.  There  we  show 
that,  under  P,  r(t)  is  (still)  strongly  consistent  and 

»  1, 

(1.45)  /t  (r  ( t)  -  r)  — >  oN(  0, 1  )/(EaQ)/2 
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where 


(1.46)  a2  •  EYQ!f)2  *  2  l  EY0(f)Y  (f)  . 

k”1 

This  CLT  corresponds  to  (1.28)  with  E  (expectation  under  P)  replacing  E.  Estimate 
2 

o  and  construct  confidence  intervals  just  as  before. 

The  stationarity  assumptions  (1.14)  and  (1.44)  can  be  significantly  relaxed.  Our 
CLT's  can  be  extended  to  Include  certain  nonstationary  processes  by  appealing  to 
Billingsley  ((1968),  Theorem  20. 2J.  Furthermore,  our  discussion  carries  over  to  the 
steady-state  estimation  problem  for 

r  -  lim  ~  ).  f(X(U  >>I  (t) 

t+m  t  k-0  *  l°'V 

(when  the  limit  exists),  where  the  CJy'e  are  an  increasing  sequence  of  random  times.  Such 
limits  are  of  interest,  for  example,  in  queues  where  lump-sum  rewards  are  payed  out  to  the 
server  at  customer  departure  epochs.  Our  arguments  go  through  provided  that  one  modifies 
the  definition  of  Yn(f)  to 


Y  (f)  -  l  f(X(0  ))I.  .(U.)  . 

n  j-0  j  Sn'Sn+1  j 

1.4  Transactions.  To  essentially  every  time  average  r  there  corresponds  a  transaction 
(custoaier)  average  s  and  conversely.  Heyman  and  Stidham  (1980)  and  Heyman  and  Sobel 
((1982),  §11.3)  establish  this  correspondence  explicitly.  Thus,  every  estimator  of  r 
yields  an  indirect  estimator  of  8  and  conversely.  In  our  framework.  Section  4  indicates 
that  on  balance  estimating  s  indirectly  is  often  better  than  estimating  it  directly. 

This  conclusion  appears  contrary  to  some  folklore. 

Readers  wishing  to  skip  to  Section  4  can  do  so  without  loss  of  continuity. 


2.  RESULTS  MID  PROOFS  I  I. 


Throughout  this  section  we  assume  that  (1.14),  (1.17),  and  (1.18)  hold  under  P. 


THEOREM  1.  Formuias  (1.23),  (1.24),  and  (1.2S)  hold  a.s. 


Proof.  First,  observe  that  (Y_(f),a  )  is  trivially  a  functional  of  X  (•)  and  thus 
— -  n  n  n 

( (Y  (f),a  i  :  n  >  0}  is  *-mixing  with  the  same  mixing  coefficients  as  the  X  's.  Since 
n  n  n 

any  4-mixing  sequence  is  ergodic  (Lamperti  (1977),  pp.  95-96),  apply  Birkhoff's  ergodic 
theorem  (Heyman  and  Sobel  (1982),  p.  366)  to  conclude  that 

(2. 1)  Y  (f )  ♦  EY  (f )  a.s. 

n  u 

(2.2)  a  ♦  Ea„  a.s. 

n  0 

proving  (1.23)  with  r  given  by  (1.25).  With  N(t)  defined  by  (1.33), 

(2.3)  SN(t)/N(t)  *  ^(t)  <  SN(t)+1/Nlt) 

from  the  definitions.  From  N(t)  t  «•  a.s.  and  (2.2),  the  extreme  terms  of  (2.3)  converge 
to  EaQ  a.s.;  as  t/H(t)  gets  squeezed,  it  converges  to  Eafl  a.s.  For  nonnegative  f, 


(2.4) 


'  N(t) 
N(t) 


W,,f) 


<  r(t) 


< 


1  (H(t>  ♦  1)  c 

t  N(t)  N(t)+1 


(f)  . 


Apply  N(t)  1  “  a.s.,  (2.1),  and  N(t)/t  ♦  1/Ea0  a.s.  to  quickly  see  that  the  extreme 

terms  of  (2.4)  converge  to  r  given  by  (1.25),  proving  (1.24)  for  f  »  0.  Split  a 
general  f  into  its  positive  and  negative  parts,  apply  (1.24)  to  each,  and  recombine.  To 
justify  the  last  step,  use  (ElY^lf)!)2  <  EY^(f)2  <  •  by  (1.18). 


THEOREM  2.  Formula  (1.26)  holds. 


Proof.  Observe  that 


,  n+1  ,  n 

<2,5)  iil  Ilx,«)(aN(s)  ,dS  *  n  J-  ‘Vtx,-)  (“k  *  *  E*ao'  a0  *  x)  a’8, 

U  K*  0 
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(E{X>A)  4  EXIa  where  is  1  or  0  depending  on  whether  or  not  A  occurs)  and  then 

mimic  the  proof  of  Theorem  1. 


THEOREM  3.  Formulas  (1.27)  and  *.28)  hold. 


Proof .  The  sequence  {Y^If)}  is  ♦-mixing  with  zero  mean.  It  satisfies  the  conditions  of 
Billingsley  (1968),  Theorem  20.1,  and  hence 

(2.6)  /n  5  (f )  — >  CN( 0,1) 

n 

as  n  ♦  •*(  for  o  ■  0,  see  Billingsley  (1968),  p.  177.  Rewrite  (2.6)  as 

(2.7)  /n  a  (r  -  r)  — >  aN(0,1)  . 

n  n 

Now  (1.27)  follows  from  (2.7  ,  (2.2),  and  a  routine  application  of  the  converging-together 
lemma  (Billingsley  (1968),  p.  25). 

For  (1.28),  let  {tn>  be  an  arbitrary  sequence  converging  to  •».  Apply  Billingsley's 
[(1968),  p.  146]  random  time  change  theorem  to  the  weak  invariance  principle  version  of 
(2.6)  to  get 

(2.8)  >rN(tk)  ?N(t  }(f)  — >  «N<0,1) 


as  k  ♦  «•.  Another  application  of  the  converging-together  lemma  then  yields 

J  -  r)  — >  ON<0.1)//^ 


using  N(t)/t  ♦  EOq 

(2. 10) 


proved  just  after  (2.3)  and  /Ntt^)  *  /tk  /N(tk )/tk  . 
^k  'r«(tk)  ”  r(tk)l  *  ’^k  YN(tk)(|f  1  ,/SN(tk)  ' 


Clearly 


2 

Combining  a  standard  Borel-Cantelli  argument  and  EYk(|f|)  <  •  proves  that 
(2.11)  Yk<|f|)//k  ♦  0  a.s. 


But  the  right-hand  side  of  (2.10)  can  be  rewritten  as 


(^(^,)2*(Ymtk)<|f|,/’^V,(H(tk>/sN(tk)) 


which,  by  (2.11),  converges  to 


(Ea„)/Z  •0»(Ean)  -  0. 


Apply  the  converging-together  lemma 
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to  (2.9)  and  (2.10)  to  get  (1.24)  for  t  going  to  »  through  the  sequence  tfc.  Because 
that  sequence  is  arbitrary,  (1.24)  holds  without  qualification  (see  Billingsley  (196B),  p. 
16). 

Now  also  assuming  (1.34)  we  have 


THEOREM  4.  Formula  (1.35)  holds  with 


(2.12 


Vo<f| 


*  1 

k«1 


(Ea0vk(f) 


+  EYgff)^) 


Proof.  The  infinite  series  for  8  converges  absolutely  by  a  remark  of  Billingsley  (1968), 
p.  177.  Turning  to  the  bias  expansion  itself,  observe  that 

(2.13)  r  -  r  -  7  (f)g(o  > 

n  n  n 

where 

(2.14)  g(x)  -  1/x  . 

Choosing  c  small  enough, 

(2.15)  sup{g"(x)  :  |x  -  EaQ |  <  e}  »  M  <  -  . 

Letting  event  Bn  -  (la^  -  Ea0l  <  e) ,  use  (1.34)  and  then  Chebyshev's  inequality  to  get 

(2.16)  |E{r  -  ri  Bc)  |  <  KP{|5  -  E«U  >  e) 

n  n  n  u 

<  KE(o  -  Eo-lVe4  . 
n  u 


By  Billingsley  ((I960),  lemma  4,  p.  172], 

(2.17)  E(an  -  Ea0)4  -  0(n"2)  . 

Next 

E{r  -  r)  *  E{r  *  ri  B  ]  *  E{r  -  ri  Bc)  . 
n  n  n  n  n 

Combining  (2. 16)— (2. 17)  gives 

(2.18)  E{r  -  r>  «  E{r  -  ri  B  >  +  0(n  2)  . 

n  n  n 

On  the  event  B„,  expand  g(an>  In  *  Taylor  series  around  BaQ  to  get 

(2.19)  E{ r  -  ri  B  }  «  S  +  T 

n  n 


where 


r 


(2.20)  s  -  «{v‘f»f57  -  -7-2  <5„  -  “VW 

0  (Ea0) 

(2.21)  T  -  «{?  <f)g*(C  Ha  -  Baa)2/2iB  } 

n  n  n  u  n 

and 

(2.22)  I5n  -  Ba0t  <  c  . 

Apply  Cauchy-Schwart*  and  than  Bllllngslay  ((1968),  Leases  3  and  4,  p.  172]  to  gat 

(2.23)  T  «  HB1 Y  (f)2)1/2(I(S  -  Ia.)4>1/2  -  0(n“3/2)  . 

n  no 

An  argument  alnllar  to  that  justifying  (2.18)  gives 

(2.24)  8  -  -*(Yn(f)Sn)/(*»0)2  ♦  0(n”2)  . 

A  alapla  calculation  shows  that 

(2.25)  B(?  (f )a  )  ■  t/n  ♦  o( 1/n) 

n  n 

with  6  given  by  (2.12).  Combining  (2. 18)— (2.25)  glvaa  (1.35),  finishing  the  proof. 

COROLLARY.  Jackknifing  worksi  (1.38)  holds. 

THBORSM  5.  Formula  (1.40)  holds. 


Proof.  By  the  converging-together  la— a,  It  suffices  to  prove  that 
( 2 * 2® )  <r2n  "  *2n>  *  0 

In  probability.  Straightforward  algebra  shows  that  the  left  side  of  (2.26)  equals 
(2.27)  ®n (Cn  *  V/<Cn  ♦  D„> 

where 


B„  •  r(0,n  -  1)  -  r(n,2n  -  1) 
C„  -  «(0,n  -  1)/n 
Dn  *  a(n,2n  -  1)/n 
b 

a(a,b)  -la,. 

j-a  3 

An  argument  He  proof  of  (1.27)  shows  that 

(2.28)  /n  Bn  — >  /2  sBIO.D/to,,  . 
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Clearly 


(2.29)  Cn  -  Dn  *  0 

Combining  (2.27)  through  (2.29)  verifier  (2.26). 


a. a. 
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3.  RgflPLTS  AND  P ROOfS:  II. 

Theorem  6.  Aaauae  that  (1.14),  (1.17),  and  (1.18)  hold  under  V.  Than,  thara  exiata  a 
probability  p  auch  thati 

,  t  . 

1)  J  J  P{X  8  e  «}da  ♦  P{X  e  •}  aa  t  ♦  • 

*  0  * 

il)  p{x  •«*•>■  P(X  e  •)  for  u  >  0 

m  ^ 

ill)  ?{X  e  •)  -  ~  !  P{X  •  8  e  *|8  >  a}da< 

BS,  0  * 

Proof .  bat  f  ba  a  bound# <3  nonnegative  function  on  £).  For  X  e  0,  aat 
fu(X)  •  f(X  •  8u)  and  put 

S1 

,U(X)  -  J  fu(x  .  eB)d.  -  ,0(X  •  eu) 


Proceeding  aa  in  the  proof  of  Theorem  1,  Birkhoff'a  argodic  thaoram  proven  that 


(3.1) 


j  S  +u 

09»<x*  V*/n  f(x  *  v*  ♦  *«<« 

K«0  k  u 


P  a. a.,  aa  n  •»  «■.  The  "squeeze"  argument  of  (2.4),  applied  to  (3.1),  can  ba  readily 
adapted  to  ahow  that 


,  u+t  1  * 

l  f  f(X  •  8  )da  ♦  —  >g„(X  •  8  )  P  a. a. 

U  BS^ 


It  followa,  from  bounded  convergence,  that 


t  . 


I  J  Bf(X  •  8a)da  ♦  tg0(X  •  6U) 


aa  t  •*  «i  in  particular,  apacialiaing  f  to  indicator  functione  yialda 
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(3.2) 


as  t  ♦  «.  Setting  u  “  0  in  (3.2)  proves  i)  and  iii).  For  ii),  observe  that  the  left- 
hand  side  of  (3.2)  is  independent  of  u. 

Throughout  the  remainder  of  this  section,  we  assume  that  (1.14),  (1.17),  and  (1.18) 
are  in  force  for  P»  and  that  P  denotes  probability  under  the  measure  P  constructed  in 
Theorem  6. 


THEOREM  7. 

(3.3)  P{r(t)  ♦  r}  -  1  . 

Proof.  Use  definitions  and  then  the  fact  that  (1.24)  holds  under  P  to  get 

(3.4)  P{r (t )  *8  ♦  ri  S.  >  s> 

B  i 

*  8+t  * 

«  p{J  f (X(u) )du/t  ♦  r»  S^  >  s|  -  P{s1  >  s}  . 

s 

Now  integrate  formula  (3.4)  with  respect  to  s  from  0  to  •».  On  the  right  we  get  ES^ , 

positive  by  our  assumptions.  On  the  left,  use  dominated  convergence  to  take  the  limit  with 

x  * 

respect  to  t  inside  the  integral  and  then  (1.13)  to  get  P[r(t)  ♦  rlES^  Cancelling 
ES1  from  both  sides  gives  (3.3). 

THEOREM  8.  Under  P  formulas  (1.45)  and  (1.46)  hold. 

Proof.  For  any  x, 
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(3.S) 


P{/t  (r(t)  -  r)  •  <  «)  S,  >  ■> 

.  a+t  . 

-  p{(1//t  j  f (X(u))  <  x|S^  >  a}{S,  >  a)du  . 

a 

By  Billingsley  ((1968).  Theorem  20.2)  the  conditional  probability  above  convergea  to  the 
appropriate  normal  probability  aa  t  ♦  «,  for  every  x  and  a.  Integrate  both  aidea  of 
(3.5)  with  respect  to  a  from  0  to  uae  (1.13)  to  aaa  that  the  left  aide  equala 

p(/t  (r(t)  -  r)  <  x}ESj,  uae  doaiinated  convergence  on  the  right  to  take  the  limit  with 
respect  to  t  inside  the  outer  integral,  and  recall  that  any  nonnegative  random  variable 

m 

2  has  expectation  J  p(z  >  t)dt.  Cancel  88.  from  both  sides  of  (3.4)  integrated  aa 
0 

above  to  get  that  p{/t  (r(t)  -  r)  <  x)  convergea  to  the  appropriate  normal  prooability. 
A  routine  calculation  verifies  (1.46). 
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4.  INDIRECT  ESTIMATION  OF  TRANSACTION  AVERAGES . 

Suppose  a  performance  measure  s  aggregates  costs  for  transactions  moving  through  the 
system.  Let  transaction  i  arrive  at  time  and  leave  at  time  Bi .  Put 

Di  “  Bi  '  Ai*  Assume  that  the  cost  associated  with  transaction  i  is  g^fD^)  and  that 
the  average  cost 

(4.1)  s  -  lim  {g,(D  >  +  .  .  +  g  (D  >)/n 

11  n  n 

exists  a.s.>  clearly 

(4.2)  sn  -  19,(0,)  +  ...  +  9n(Dn)]/n 

consistently  estimates  s.  Our  discussion  carries  over  to  more  general  transaction 
averages.  He  use  s  for  concreteness. 

In  the  notation  of  Section  1,  state  x^  includes  the  value  of  for  every 

transaction  i  in  the  system  at  time  T^  and 

(4.3)  f(X(t))  -  l  g*(t  -  A1) 

where  the  sum  is  over  those  transactions  i  for  which  A^  <  t  <  .  Here  gj  is  the 

right  derivative  of  gA  and  by  definition  g^(d!  -  0  for  d  <  0  and  for  d  >  D^,  This 

setup  allows  rather  general  g^’s,  though  not  indicators. 

Now  call  the  arrival  rate  X.  Assuming  all  limits  exist,  r  »  Xs  by  Heyman  and 

Stidham  (1980).  Usually  X  la  simply  the  reciprocal  of  the  expected  arrival  spacing.  In 

more  general  cases  such  as  batch  arrivals  with  batch  size  and  spacing  dependent,  the  theory 

developed  earlier  in  this  paper  shows  how  to  estimate  X. 

He  can  choose  whether  to  observe  discretized  observations  y  (f)  as  detailed  earlier 

n 

or  transaction  observations  g^  (D^ ) .  Our  choice  is  based  on  four  criteria: 

(1)  Ease  of  data  collection.  Formula  (4.3)  indicates  that  gathering  discretized 

observations  at  state-change  epochs  is  more  work,  but  probably  not  much  more  because  the 

evolution  of  the  system  has  to  be  simulated  in  any  case.  Gathering  discretized 

observations  f(X(t))  at  (reasonably  short)  equally- spaced  intervals  that  do  not 

necessarily  correspond  to  state-change  epochs  may  be  much  more  work. 

(ii)  True  variance.  Some  believe  that  Xs  always  has  smaller  variance  than  r. 
-  n  ® 

where  n  is  the  number  of  transactions  completely  processed  up  to  event  epoch  m.  To  some 
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extent  the  results  of  Cerson  end  Law  (1980)  for  special  esses  support  this  folklore,  but 
Cooper's  ((1981),  pp-  293-296]  example  shows  that  It  is  untrustworthy. 

(Ill)  Variance  estlaatlon.  Folklore  has  It  that  the  sequence  of  transaction  observations 

{g^ (D^) }  is  covariance  stationary  if  and  only  if  the  sequence  of  discretised  observations 

(Y  (f)}  is  covariance  stationary,  thus,  loosely  speaking,  any  variance-estimation 
n 

technique  (e.g. ,  via  batch  neans,  spectral  analysis,  autoregressive  representations, 
functional  limit  theorems)  applies  to  both  or  to  neither.  In  practice,  however,  one  has  to 
ask  at  what  sample  sixes  asymptotic  results  reasonably  apply)  the  faster  covariance  falls 
off  with  increasing  lag,  the  better.  Generally,  for  discretized  observations,  covariance 
does  drop  fairly  quickly.  By  contrast,  especially  in  systems  that  are  not  first-in,  first- 
out  (FIFO),  a  series  of  transaction  observations  scrambles  past,  present,  and  future  and 
cuts  connections  between  time  spacing  and  index  spacing.  Thus,  covariance  between  widely- 
spaced  observations  may  well  be  significant.  This  messy  covariance  structure  is  hard  to 
handle,  at  least  with  practical  samaple  sizes. 

(lv)  Bias.  The  scrambling  affect  mentioned  in  point  (iii)  above  may  make  it  harder  to 

detect  (and  hence  attenuate)  initialization  bias  for  transaction  observations. 

Incompletely  processed  transactions  cause  termination  bias  for  transaction  observations, 

exacerbated  by  the  "inspection  paradox"  discussed  earlier.  For  discretized  observations, 

r  is  generally  biased  but  our  jackknife  using  r  makes  this  bias  nearly  negligible. 

**  n 

We  conclude  that,  at  least  for  non-FIFO  systems,  indirect  estimation  via  discretized 
observations  has  more  advantages  than  disadvantages  relative  to  direct  estisuttion  via 
tranaaction  observations.  The  results  developed  in  the  preceding  sections  make  statistical 
analysis  of  the  former  possible  and  practical. 

Asymptotically,  s  -  r/X.  Estimate  s  by  plugging  in  consistent  estimators  for  r 
and  X.  If,  as  is  usual,  X  is  known,  then  divide  confidence  limits  for  r  by  X  to  get 
confidence  limits  for  s.  If  X  has  to  be  estimated,  a  straightforward  ratio-estimator 
analog  of  the  methods  discussed  in  this  paper  must  be  used  to  obtain  an  asymptotically 
exact  confidence  interval. 


-19- 


REFERENCES 


1.  P.  BILLINGSLEY  (1968),  Convergence  of  Probability  Measures,  Wiley,  New  York. 

2.  P.  BRATLEY,  B.  L.  FOX,  and  L.  E.  SCHRAGE  (1983),  A  Guide  to  Simulation, 
Springer-Verlag,  New  York. 

3.  J.  S.  CARSON  and  A.  M.  LAW  (1980),  Conservation  equations  and  variance  reduction  in 
queuing  simulations,  Oper.  Res.,  28,  pp.  535-546. 

4.  R.  B.  COOPER  (1981),  Introduction  to  Queueing  Theory,  2nd  edn..  North  Holland,  New 
York. 

5.  P.  FRANKEN,  D.  KONIG,  0.  ARNDT,  and  V.  SCHMIDT  (1982),  Queues  and  Point  Processes, 
Wiley,  New  York. 

6.  P.  W.  GLYNN  and  D.  L.  IGLEHART  (1981),  Simulation  output  analysis  for  general  state 
space  Markov  chains,  Proc.  of  the  ORSA-TIMS  Special  Interest  Meeting  on  Applied 
Probability  -  Computer  Science,  The  Interface. 

7.  D.  P.  HEYMAN  and  M.  J.  SOBEL  (1982),  Stochastic  Models  in  Operations  Research,  Vol.  1, 
McGraw-Hill,  New  York. 

8.  D.  P.  HEYMAN  and  S.  STIDHAM,  Jr.  (1980),  The  relation  between  customer  and  time 
averages  In  queues,  Oper.  Res.,  28,  pp,  983-994. 

9.  L.  JOW  (1982),  An  autoregressive  method  for  simulation  output  analysis.  Technical 
Report  21,  Dept,  of  Operations  Research,  Stanford  University,  Stanford,  California. 

10.  J.  LAMPERTI  (1977),  Stochastic  Processes.  Springer-Verlag,  New  York. 

11.  D.  R.  MILLER  (1972),  Existence  of  limits  in  regenerative  processes,  Ann.  Math.  Stat.  , 
43,  pp.  1275-1282. 

12.  A.  STRELLER  (1980),  A  generalization  of  cumulative  processes,  Elektron. 
Informationsverarb.  Kybernetik,  16,  pp.  449-460. 

13.  R.  WOLFF  (1982),  Poisson  arrivals  see  time  averages,  Oper.  Res.,  30,  pp.  223-231. 


BLF/PWG/ed 


SECURITY  CLASSIFICATION  OF  THIS  RAGE  (Whn  Dmtm  Rntf4) 


s  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

S.  RECIPIENT'S  CATALOG  NUMBER 

/S' 

«.  TITLE  (end  Subtitle) 

ESTIMATING  TIME  AVERAGES  VIA  RANDOMLY- SPACED 
OBSERVATIONS 

S.  TYPE  OP  REPORT  B  PERIOD  COVERED 

Summary  Report  -  no  specific 
reporting  period 

S.  PERPORMING  ORG.  REPORT  NUMBER 

7.  AUTHORfa) 

Bennett  L.  Fox  and  Peter  W.  Glynn 

B.  CONTRACT  OR  GRANT  NUMBERfA) 

DAAG29-80-C-0041 

#.  PERPORMINO  ORGANIZATION  NAME  AND  AOORESS 

Mathematics  Research  Center,  University  of 

610  Walnut  Street  Wisconsin 

Madison.  Wisconsin  53706 

10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Work  Unit  Number  5  - 
Optimization  and 

Large  Scale  Systems 

>1.  CONTROLLING  OPPICE  NAME  AND  ADDRESS 

U.  S.  Army  Research  Office 

P.O.  Box  12211 

Research  Triangle  Park,  North  Carolina  27709 

12.  REPORT  DATE 

January  1984 

IS.  NUMBER  OP  PAGES 

20 

14.  MONITORING  agency  name  B  AODRESyif  dl  He  tent  Bom  Controlling  Olllce) 

IS.  SECURITY  CLASS,  (ol  title  report) 

UNCLASSIFIED 

IB*.  DECLASSIFY  ATI  ON/ DOWN  GRADING 
SCHEDULE 

IS.  DISTRIBUTION  STATEMENT  (ol  M*  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (ol  Ike  mbatrmel  entered  In  Block  30,  II  dWdront  Boot  Report) 


IB.  SUPPLEMENTARY  NOTES 


IB.  KEY  WORDS  (Continue  on  rororto  e  Ido  II  nocoooooy  and  Identity  by  block  number) 


Simulation,  point  processes 


r20TAIISTRACT(Coniinueonrerereeeldellneceeeaiyiiidldentifybybiocklnimberj 

To  estimate  continuous-time  averages  via  randomly-spaced  observations  of 
discrete-event  systems,  we  develop  a  point-process  framework  and  use  it  to 
generalize  both  regenerative  and  stationary-process  oriented  simulation 
methodologies.  We  give  consistent  estimators,  central  limit  theorems,  and  an 
effective  bias-reducing  jackknife.  The  impact  on  indirect  estimation  of 
transaction  (customer)  averages  is  discussed. 
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